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ABSTRACT 

A stalled spherical accretion shock, such as that arising in core-collapse supernovae, is unstable to non- 
spherical perturbations. In three dimensions, this Standing Accretion Shock Instability (SASI) can develop 
spiral modes that spin-up the protoneutron star. Here we study these non-axisymmetric modes by combin- 
ing linear stability analysis and three-dimensional, time-dependent hydrodynamic simulations with Zeus-MP, 
focusing on characterizing their spatial structure and angular momentum content. We do not impose any ro- 
tation on the background accretion flow, and use simplified microphysics with no neutrino heating or nuclear 
dissociation. Spiral modes are examined in isolation by choosing flow parameters such that only the funda- 
mental mode is unstable for a given polar index I, leading to good agreement with linear theory. We find that 
any superposition of sloshing modes with non-zero relative phases survives in the nonlinear regime and leads 
to angular momentum redistribution. It follows that the range of perturbations required to obtain spin-up is 
broader than that needed to obtain the limiting case of a phase shift of ir/2. The bulk of the angular momentum 
redistribution occurs during a phase of exponential growth, and arises from internal torques that are second 
order in the perturbation amplitude. This redistribution gives rise to at least two counter rotating regions, with 
the maximum angular momentum of a given sign approaching a significant fraction of the mass accretion rate 
times the shock radius squared (Mr s 2 hl , k ~ 10 47 g cm 2 s , spin period ~ 60 ms). Nonlinear mode coupling at 
saturation causes the angular momentum to fluctuate in all directions with much smaller amplitudes. 
Subject headings: hydrodynamics — instabilities — shock waves — supernovae: general — stars: rotation — 
pulsars: general 



1. INTRODUCTION 

The explosion mechanism of massive stars is not well un- 
derstood at present. Observational and theoretical evidence 
gathered thus far suggests that this mechanism is intrinsi- 
cally asy mmetric for all stars that form iron cores (> \2M p,\ 
see, e.g.. IWang & Wheeler! 120081 iBurrows et aD 12007b! and 
iJanka et all 12007 for recent reviews). Stars that end up with 
O-Ne-Mg cores are currently found to explode in spherical 
symmetry via the neutrino mechanism (jKitaura et al. 2006, 
IBurrows eTaT]|2007al IJanka et al.ll2008l) . 

One piece of phenomenology that a successful theory of 
core-collapse supernovae has to explain is the distribution of 
pulsar spins at birth. Population synthesis studies of radio 
pulsars generally assume normal or log- normal distributions 
with m ean values ranging from a few ms dArzoumanian et alj 
120021) to a few ~ 100 ms ( Faucher-Giguere & Kaspill2006l) in 
order to reproduce observations. This large range is due to 
the different assumptions made abo ut input physics, such as 
the sh ape of the radio beam (e.g., Faucher-Gig uere & Kaspil 
2006). Independent constraints that combine Chandra, XMM, 
and Swift observations of historic supernovae with an em- 
pirical correlation between X-ray luminosity and spin down 
power tend to rule o ut a large pulsar population with ini- 
tial periods < 40 ms dPerna et alJl2008l) . On the other hand, 
current stellar evolution calculations that account for mag- 
netic torques predict, using conservation of angular momen- 
tum and an assump tion for the mas s cut, neutron star birth 
periods < 10 ms (Hegeretal. 2005). However, these mod- 
els suffer from large uncertainties in the input physics, which 
could also lead in princ iple to very slowly rotating pulsars 
(Sprui t~& Phinneylll998h . Axisymmetric core-collapse calcu- 
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lations indicate that there is an almost linear mapping between 
the rotation rate of the iron core and that of the resulting pro- 
toneutron star, with no robust braking mechanism in sight to 
bring the implied fast neutro n star spins to values more in 
agreement with observations dOtt et al.ll2006l) . 

If most presupernova cores turn out to rotate slowly, then 
there is an alternative mechanism to generate periods > 
50 ms, which arises from instabilities in the supernova shock 
itself (Blondin & Mezzacappa] 120071) . Axisymmetric core- 
collapse supernova simulations that include neutrino trans- 
port and microphysics to several degrees of approxima- 
tion find that the stalled postbounce shock undergoes large 
scale oscillator y motions that break the spherical symmetry 
of the system (IBurrows ^ et all 1 19951: Uanka & Mueller! 1 19961: 
Mezzacappa et all 119981: IScheck et all 120061: lOhnishi et al] 
20061: iBuras et al] I2006bllat IBurrows et al] 120061 I2007d: 
Scheck etaljl20 08: Murp hy & Burrowsll2008l:IOtt et alJl200a 
Marek & Jankal 120091: ISuwa et all 12009^ When neutrino 
driven convection is suppressed, the instability of the shock 
persists in the form of an overstable sloshing cycle, with the 
most unstable mode s having Legendre indices I = 1 and 2 
dBlondin et all 120031: iBlondin & Mezzacappa! 120061) . In the 
limit of short wavelength, this so-called Standing, Spheri- 
cal, or Stationary Accretion Shock Instability (SASI) is driven 
by an interplay betwe en advected and acoustic perturbations 
(Fogli zzo et al.l f2007). There is no conclusive proof yet on 
the mechanism behind long wavelength modes, although con- 
siderations of the timescales (Fernandez & Thompson 2009b) 
and t he case of planar geometry (Foglizzo 2009j ISato et all 
2009) point to the advective-acoustic cycle as a relevant com- 
ponent. The SASI is we akened when r otation becomes dy- 
namically important (e.g. JOtt et al.ll2008l) . 

Three-dimensional simulations of the SASI have found that 
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a spiral type of mode can develop, causing the flow to di- 
vide itself into at least two counter-r otating streams, leading 
to ang u lar momentum redistribu tion (Blondin & Mezzacappa 
2007). iBlondin & Mezzacappa! (120071) find that this process 
can spin-up a canonical neutron star to a period of ~ 50 ms 
after ~ 0.1M© of material is accreted. Aside from the fact 
that this spin period is comparable to that obtained by pop- 
ulation synthesis calculations, the most interesting results of 
IBlondin & Mezzacappa (2007) are that (i) no progenitor ro- 
tation is required for this mechanism to operate and, (ii) for 
weak progenitor rotation, this instability dominates the spin- 
up of the protoneutron star, imparting it along a different axis. 

One of the key questions that remains to be answered 
is how easy it is to excite these sp iral modes. In con- 
trast to IBlondin & Mezzacappa! (120071) and IBlondin & Shaw l 
( 20071). the fi ndings of llwakami et al.l d2008al ). Ilwakami et al.l 
d2009al) . and llwakami et alj (|2009b) suggest that with no ro- 
tation in the upstream flow, th ese modes are difficult to excite. 
Yamasaki & Foglizzo (2008) find that, indeed, for a cylindri- 
cal accretion flow with rotation, prograde spiral modes have 
a higher growth rate than retrograde ones. They attribute this 
to a Doppler shift of the mode frequency induced by rotation. 
Part of the difficulty in comparing numerical results on this in- 
stability is that they have been obtained using different micro- 
physics, numerical methods, and coordinate systems, and that 
the three-dimensional structure of these linea r eigenmodes re- 
mains as of yet unexplored. The work of Blondin & Shawl 
(120071) examined spiral S ASI modes with time-dependent hy- 
drodynamical simulations restricted to a polar wedge around 
the equatorial plane. They found that sloshing and spiral 
modes are closely related, because the pressure perturbation 
rotates with the same frequency as a sloshing mode of the 
same Legendre index, and because by combining two counter- 
rotating spira ls, the flow field of a slo shing mode is recovered. 
The work of Iwaka rm et al.l (I2009bl) also identified sloshing 
modes as the superposition of two counter-rotating spirals. 

The aim of this paper is to better understand spiral modes 
in the linear and nonlinear phase, particularly their three- 
dimensional spatial structure and angular momentum content. 

The structure of axisymmetric sloshing eigenmodes has 
bee n obtained previousl y through linear stability analy- 
sis ((Foglizzo et al. 2007), and verified to high precision 
with t ime-dependent axisymmetric h ydrodynamic simula- 
tions (Fernandez & Thompson 2009b). Starting from the 
findings of Blondi n & Shawl (120071) . we show that spiral 
modes are most easily understood as sloshing modes out of 
phase. Their spatial s tructure is buil t using solutions to the 
differential system of iFoglizzo et al.l (120071) . and their evo- 
lution is compared with results of time-dependent h ydrody- 
namic simulations using Zeus-MP (lHaves et al.ll2006h . In our 
calculations we do not add rotation to the flow anywhere, and 
employ an ideal gas equation of state, parametric cooling, and 
point-mass gravity, with no heating or nuclear dissociation be- 
low the shock. 

One of our main findings is that it is relatively simple to 
create spiral modes out of sloshing modes by changing their 
numbers, relative phases, and amplitudes, covering a broader 
parameter space than the limiting cases where the amplitudes 
are equal and the relative phase is tt/2. These modes survive 
to large amplitudes, resulting in a significant angular momen- 
tum redistribution below the shock. We set aside for now the 
question of whether this coherent superposition of sloshing 
modes is able to develop in a more realistic core-collapse con- 
text, where three-dimensional, nonlinear turbulent convection 



may li kely act as a forcing agent (Fernand ez & Tho mpson 
2009a). We revisit this issue in the discussion section. 

We also find that the bulk of the angular momentum redis- 
tribution generated by a spiral mode occurs during the phase 
of exponential growth. This spin-up of the flow arises from in- 
ternal torques that are second order in the perturbation ampli- 
tude, and results in angular momenta of a characteristic mag- 
nitude at saturation. Non-linear mode coupling causes the an- 
gular momentum flux to become stochastic, fluctuating in all 
directions with an amplitude much smaller than that achieved 
during exponential growth. 

The structure of this paper is the following. In Section 2 we 
describe the physical model used, and summarize the most 
important elements that enter the linear stability calculation 
and time-dependent simulations, with details deferred to the 
Appendices. Section 3 describes how spiral modes can be 
constructed by superposing sloshing modes out of phase, and 
explores some of their features. Section 4 contains the results 
of time-dependent simulations of these modes, including both 
linear and nonlinear development. We conclude by summa- 
rizing our findings and discussing their implications for more 
realistic core-collapse models and neutron star spins. 

2. METHODS 
2.1. Physical Model 

The physical system consists of a steady-state, standing 
spherical accretion shock at a radial distance r s o from the ori- 
gin. Below the shock, the material settles subsonically onto 
a star of radius r* and mass M centered at the origin of the 
coordinate system. This settling is mediated by a cooling rate 
per unit volume that mimics neutrino emission due to electron 
captu re, Jzf oc p 3 / 2 p, with p the pressure and p the density 
(e.g.. IBlondin & Mezzacappall2006l IFernandez & Tho mpson 
2009b). In all the cases studied here, cooling is relevant 
only in a narrow region - of order a pressure scale height 
- outside the accretor. Upstream of the shock, the fluid is 
supersonic, with incident Mach number M. \ = 5 at r = r s o, 
and has a vanishing energy flux. Upstream and downstream 
solutions are connected by the Rankine-Hugoniot jump 
conditions. The equation of state is that of an ideal gas of 
adiabatic index 7 = 4/3, and the self gravity of the flow 
is neglected. T he eq uations describing this steady flow 
are presented in 3A.1I T he model was first developed by 
iHouck & Cheval ier (1992) to study a ccretion onto c ompac t 
objects, and subsequen tl y used by Blon din et al.l d2003|) . 
I Blondin & Mezzacappa] (12006b. IFoglizzo et al.l ^20071) . 
Blondin & Mezzacappa! d2007l). I Blondin & Shawl d2007l) . and 
IFernandez & Thompson ( 2009bJ) as a minimal approximation 
to the postboun ce stalled shock. 

In contrast to Fernandez & Thompson (2009b), we do not 
include here the effects of nuclear dissociation at the shock. 
This is a significant energy sink in the system, changing the 
Mach number and density below the shock for fixed upstream 
conditions, and hence affecting the linear modes of the SASI. 
Nevertheless, as our goal here is to characterize the basic be- 
havior of the instability in three dimensions, we aim at keep- 
ing the calculations as simple as possible. 

Neutrino heating is also neglected, in ord er to suppress 
convec tion in the flow, as originally done by Blondin et al. 
d2003l) . Since most of the flow below the shock is adiabatic, 
the entropy profile is flat and hence marginally unstable to 
overturns triggered by the SASI itself. 

We adopt a system of units normalized to the shock ra- 
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dius r s o, the free-fall velocity at this radius vgo = y/2GM /r 5 o, 
and t he upstream density p\ (Fernan dez & Tho mpson 
2009b). Numerical values relevant for the stalled shock 
phase of core collapse supernovae are r s o ~ 150 km, 

Vff (r s0 ) ~ 4.8 x 10 9 Af};f (rso/150 km)" 1 / 2 km s" 1 , f ff = 

rso/vff - 3.1Af7_3 /2 (r s0 /150 km) 3 / 2 ms, and pi - 4.4 x 

10 7 .Mo.3^i3 ( r ,rf)/150 km)" 3 / 2 g cm" 3 (assuming a strong 
shock), where the mass accretion rate has been normal- 
ized to M = O.3MO.3M0 s" 1 . We will express quanti- 
ties involving angular momentum in terms of Mr s o ~ A-tt x 
10 46 (r s0 /150 km) 2 M .3 g cm 2 s" 1 . 

2.2. Linear Stability Analysis 

To obtain linear SASI eigenmodes and eigenfrequen- 
cies, we solve the differential system of iFoglizzo et al.l 
d2007l) . the results of which agree with time-dependent 
axisymmetric hydrodynamic s imulations to high precision 
(Fernandez & Thompson 2009b). The relevant equations are 
described in detail in Appendix |Aj here we summarize the 
main elements of the calculation. The eigenvalue problem 
is formulated in terms of perturbation variables that are con- 
served when advected. These are the perturbed mass flux 
h = Sp/ p + 5v r /v r , Bernoulli flux / = v r 8v r + Sc 2 /(j— 1), en- 
tropy 5S = (7-1)" 1 ISp/p—^Sp/p], and an entropy-vortex 
combination 5K = r 2 \ ■ (V x 5w) + £(£ + l)(p/p)5S (see also 
eq. 1A33I ). with c 2 = 7 p / p, v the velocity, and w the vortic- 
ity. Perturbations are decomposed in Fourier modes in time 
and spherical harmonics in the angular direction, having the 
general form 

6q = 6q(r)e-^Y?(e,4,), (1) 

where 5q(r) is the complex amplitude, Yf is a spherical har- 
monic, and the complex frequency satisfies 10 = w osc + /cj grow . 
This choice results in all thermodynamic variables along with 
6v r being p roportional to F„ m (cor responding to spheroidal 
modes, e.g. lAndersson et al.ll2.Q01b . The transverse compo- 
nents of the velocity involve angular derivatives of the F™, 
and have the same radial amplitude (Appendix lAl. 



5v8 = 6vn(r)e- iut -^Yt m (9, ( f>) 

Sv^ = 5v n (r)e-^' J- -jL Yf(0, 0). 
sinfl ocp 
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The system of ordinary differential equations that deter- 
mines the complex amplitudes is independent of the azimuthal 
number m of the mode. Boundary conditions at the shock 
are expressed in terms proportional to the shock displace- 
ment A£, or the shock velocity Av = — UoA£. By imposing 
Sv r (r*) = at the inner boundary, a complex eigenvalue ui is 
obtained for a given set of flow parameters. For each t, a 
discrete set of overtones is obtained, wh ich are related to the 
number of nodes in the radial d irection (Fog lizzo et alj|2007t 
Ferna ndez & Thompson! l2009bl) . The resulting eigenmodes 
have both real and imaginary eigenfrequencies. 

For fixed 7, and functional form of the cooling func- 
tion, the eigenfrequencies of the SASI depend only on the rel- 
ative size of the shock and accreting star, quantified by the ra- 
tio r*/r S Q. Time-dependent studies of individual SASI modes 
are possible by choosing model parameters such that only the 
fundamental mode is unstable, for a given £. Based on the fact 
that the relevant modes in more realistic core-collapse simula- 



tions are £ = 1, and 2, we adopt two configurations for single - 
mod e studies: r„/V s o = 0. 5 for 1=1, and r„ j r s p = 0.6 for £ = 2 
(see IFoglizzo et al.ll2007l and iFernandez & Thompsonl 2009b 
for more extended parameter studies of the SASI eigenfre- 
quencies). In addition, we explore a configuration with a 
larger postshock cavity, r*/r s o = 0.2, to more closely resemble 
actual stalled supernova shocks. In this case unstable over- 
tones of £ = 1 and £ = 2 are present. 

To facilitate visualization, and unless otherwise noted, we 
employ the real representation of spherical harmonics. For 
£=l, one has 



(4) 
(5) 
(6) 
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This basis is entirely equivale nt to the one involvin g complex 
exponentials in azimuth (e.g., Arfken & Weber 2005), but has 
a straightforward interpretation: equations (01, (0, and (O are 
dipoles in the z-, x-, and y-directions, respectively. For £ = 2, 
we have 
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(7) 
(8) 
(9) 
(10) 
(11) 



corresponding to the usual z-symmetric quadrupole (F 2 °X an d 
a series of 4-striped "beach balls" with alternating polarity and 

symmetry axis along i (l^X $ (Y^ 1 ), and z (Y^ 2 ). 

Whenever necessary, we will denote the traditional complex 
spherical harmonics as T™. 

2.3. Time Dependent Numerical Simulations 

We perform time-dependent hydrodynamic calculations to 
verify the evolution of linear SASI eigenmodes in three di- 
mensions, covering the linear and nonlinear phases. To 
this end, we empl oy the publicly available code Zeus-MP 
dHaves et al.l 120061) . which solves the Euler equations using 
a finite difference algorithm that includes artificial viscosity 
for the treatment of shocks. The default version of the code 
supports an ideal gas equation of state and point mass gravity. 
We have extended it to account for optically thin cooling and a 
tensor artificial viscosity, for a better shock treatment in curvi- 
linear coordinates (e . g.. | Stqne & Normanl l l992l lHaves et"al] 
2006, Iwa kami et al.l I2008ah . Issues associated with imple- 
menting the latter are discussed in Appendix 151 

The initial conditions are set by the steady-state accretion 
flow used in linear stability calculations. To prevent runaway 
cooling at the base of the fl ow, we impose a gaussian cutof f 
in the cooling with entropy (Ferna ndez & Th ompson 2009b). 
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The normalization of the cooling function is adjusted so that 
the Mach number at the surface of the accreting star is M ~ 
1(T 2 . 

Spherical polar coordinates (r,9,<fi) are used , covering the 
whole sphere minus a cone of half-opening angle 5 degrees 
around the polar axis. This prescription does not signifi- 
cantly alter the flow dynamics, and ameliorates the severe 
Courant-Friedrichs-Lewy restriction around the polar axis (H- 
Th. Janka, private communication). We present convergence 
studies in Appendix iBl showing that only ~ 10% differences 
relative to using the full sphere are introduced with this pre- 
scription. Cells are uniformly spaced in the polar and az- 
imuthal directions, while ratioed in radius. The choice of grid 
spacing is key to a achieve a stable background solution and 
to adequately capture linear growth rates, while minimizing 
the computational cost. The radial size of the cells at r = r*, 
Ar m i n , determines how stable the unperturbed shock remains, 
because near hydrostatic equilibrium needs to be maintained. 
At the shock, the angular size determines how well growth 
rates are captured. Based on numerical experiments, we have 
found that Ar m ; n < 10~ 3 r s o is required to obtain a stable shock 
within 500 dynamical times, and that 36 cells in the polar di- 
rection is the minimum needed to obtain a clean shock os- 
cillation (w osc within < 10% and w grow within ~ 20% from 
the linear stability value; convergence tests are presented in 
Appendix IBl. We choose the ratio of radial spacing so that 
Ar = r s oA9 at the shock radius r s o, and A<fi = Ad, resulting in 
a total of 56 radial {£ = 1), 48 polar, and 96 azimuthal cells. 
For i = 2, we use 28 ratioed cells in radius from r» to r s o, and 
then use a uniform radial spacing until the outer boundary. 
This choice of grid spacing is made because the tensor artifi- 
cial viscosity be haves best when applied over the longest cell 
dimension (e.g., IStone & Normanlll992l) . and hence choos- 
ing the three cell dimensions as equal as possible results in 
nearly isotropic dissipation except near the polar axifl In ad- 
dition, we perform one simulation at double the resolution 
in all dimensions (e.g., 128x96x192) to test the robustness 
of our results. This resolution approaches values comparable 
to existing two-dimensional simulations, and can better cap- 
ture the parasiti c instabilities tha t mediate the saturation of 
the instability (Guil et et al.ll2010t) . In Table Q] we summarize 
the modes evolved. With the exception of the high-resolution 
model (R5_L1P2_HR), all simulations were carried out until 
at least 200fffo, well into the nonlinear phase. 

The boundary conditions are reflecting at the inner radial 
boundary (r = r*), periodic in the <\> direction, reflecting at the 
polar boundaries (6 = {0,7r} minus the cutout), and set equal 
to the upstream flow at the outer radial boundary (r = 4r s o). 
The use of tensor artificial viscosity prevents the appearance 
of a nume rical instability at the poles of the shock, as de- 
scribed byjlwakami et al. (2008b). This carbuncle instability 
dOuirklll994l)" arises within a few dynamical times if the stan- 
dard IVonNeum arin & Richtmveri (Q~950) prescription is used. 
However, even when including the tensor viscosity, we have 
found that a different numerical instability develops around 
the polar axis when reflecting boundary conditions and the 
full sphere are used. This is due to the fact that a wake is gen- 
erated by the reflecting axis when the sloshing modes reach a 
large amplitude, resulting in a different radial velocity profile 
at opposite sides of the polar axis. The azimuthal velocity then 
develops a sawtooth instability around the axis. The cutout 

2 The convergence of meridional grid lines inevitably results in a much 
smaller cell size in the azimuthal direction next to the shock. 



TABLE 1 
Models Evolved. 



Model a 


r*/r s0 


Perturbation b 


Resolution c 


—C T 1 T ~D(r-\ 


U.J 


Shell, I = 1 


56x48 


1 J_L 1 _rlK^C ) 


U.J 


Shell, i = 1 


112x96 


R2 Lllffc) 


0.2 


Shells, $ = 7r/2 


56 x 48 x 96 


R2_Lllh(c) 


0.2 


Shells, $ = tt/2 


56 x 48 x 96 


R5_Lllx 


0.5 


Shell, Fj 1 only 


56 x 48 x 96 


R5 LllP2(c) 


0.5 


Shells, $ = 71-/2 


56 x 48 x 96 


R5 L11P4 


0.5 


Shells, $ = 71-/4 


56 x 48 x 96 


R5 L11P8 


0.5 


Shells, $ = 7r/8 


56 x 48 x 96 


R5 Lll HR 


0.5 


Shells, $ = 7r/2 


112x96x 192 


R5_RAN 


0.5 


Random Pressure 1% 56 X 48 X 96 


R6 L21P2(c) 


0.6 


Shells, $ = 7r/2 


(28 + 44) x 48x96 


R6 L21P4 


0.6 


Shells, $ = 7r/4 


(28+44) x 48x96 


R6 L21P8 


0.6 


Shells, $ = 7r/8 


(28 + 44) x 48x96 


R6_L22P2(c) 


0.6 


Shells, $ = 7r/2 


(28+44) x 48x96 


R6 L22P4 


0.6 


Shells, $ = 7r/4 


(28+44) x 48x96 


R6 L22P8 


0.6 


Shells, $ = 7r/8 


(28+44) x 48x96 


a In three-dimensional runs, we denote by L£lml the case where both 
m modes are excited with the quoted phase difference $. The symbol 



(c) indicates that, in addition to the 5 degree cutout around the polar 
axis, a run with the full range of polar angles was performed in the linear 
phase only. In this case, the letter c is appended to the model name to 
denote the version with cutout. The letters / and It mean that either the 
fundamental, or first overtone was excited, respectively. 
b See eq. U2t for the definition of $. 

c Number of cells in the r, 9, and <j> directions, respectively. For the 
r*/r s o = 0.6 runs, we use a ratioed grid for r < r s o, and thereafter keep 
the radial cell spacing constant. 

improves (but does not completely suppress) the development 
of this instability, to the point where it does not significantly 
alter the flow dynamics. 

3. SUPERPOSITION OF LINEAR EIGENMODES 

Understanding of spiral SASI modes becomes easier 
when the real representation of spherical harmonics is used 
(eqns. EJ-flTTJ). For I = 1, they correspond to the familiar 
sloshing modes, axisymmetric relative to the three cartesian 
coordinate axes. When viewed this way, the fact that the 
eigenfrequencies for a given i are independent of azimuthal 
number m is natural, as any of these elementary modes is 
equally likely to arise in a spherically symmetric accretion 
flow (see, e.g.. lBinnev & Tremaindl2008l for a more rigorous 
argument). When excited in phase, any linear combination of 
these sloshing modes will result in a new sloshing mode that 
is symmetric around some axis. 

Spiral modes arise whenever two or more of these "ba- 
sis" sloshing modes are excited out of phase. The familiar 
£ = 1, m = ±1 spirals corresponds to ± /Ff 1 , or x- and y- 
symmetric sloshing modes that are ±7r/2 out of phase with 
each other, and that have the same amplitude, akin to two 
transverse and out-of-phase harmonic oscillators describing 
circular motion. Aside from the exponential growth in am- 
plitude, this gives rise to a static spiral pattern that rotates 
with the angular frequency of the oscillators. T his type of 
spiral motion was suggested by iLin & Shul (1 1964b to explain 
the morphology of spiral arms in gal axies, and is also well- 
known in stellar pulsation theory (e.g. JUnno et al.lfl989l) . 

Figure Q] shows, side-by-side, the real part of the entropy, 
entropy-vortex, pressure, and azimuthal velocity perturba- 
tion in the equatorial plane (z = 0), with the total velocity 
field superimposed. Shown are the fundamental modes Ff 1 , 
exp(/7r/2) Ff 1 , and the sum of them, for a shock displace- 
ment amplitude A£ = 0.25r s o (an animated version of the 
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FIG. 1 . — Real part of perturbed quantities from linear stability analysis, showing how transverse sloshing modes out of phase combine to create a t = 1, m = 1 
spiral mode at the equatorial plane (z = 0). The top, middle, and bottom rows show the mode at t = 0, the Y^ 1 mode displaced 7r/2 in phase, and their 
sum, respectively (real spherical harmonics are given by eqns. d-GD). Fr om left to right, columns show the entropy, entropy-vortex, pressure, and azimuthal 
velocity perturbation for a shock displacement amplitude A£ = 0.25r s o ( j|2.2l and AppendixfA)- White arrows show the perturbed velocity profile. The size of 
the postshock cavity is r*/r s o = 0.5, and only the fundamental mode is unstable. Note that only the region below the unperturbed shock position is shown. The 
electronic version of the article contains an animated version of this figure (real frequency only). 
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FIG. 2. — Effect of changing the relative phase (left) and amplitude (right) 
relative to the I = 1, m = 1 spiral in Figure Q]. Shown is the perturbed az- 
imuthal velocity from the linear eigenmodes at z = 0. The left panel has a 
phase difference <I> = 7r/8 with equal amplitudes, while the right panel de- 
creases the amplitude of Fj -1 by 1/2, with a phase difference of n/2. The 
electronic version of the article contains an animated version of this figure 
(real frequency only). 



figure is available in the online version of the article). In 
sloshing modes, the entropy eigenmodes consist of annular 
regions of alternating polarity that extend over half a circle 
in azimuth for £ = 1, and which are advected with the back- 
ground flow. The perturbed fluid moves towards regions of 
low-entropy, which generally coincide with underpressures. 
The entropy-vortex perturbation (5K) has opposite polarity 
than SS. When combined tt/2 out-of-phase, a spiral pattern is 
created in SS, SK, and Sv^. The pressure, however, has a dif- 
ferent behavior. In a sloshing mode, it propagates outwards, 
as can be seen in the animated version of Figure Q] The lat- 
ter suggests that pressure perturbations arise in response to 
negative entropy perturbations, which affect the cooling and 
thus hydrostatic equilibrium (a s in the case of £ = modes, 
iFernandez & Thompson! l2009bh . Th e source terms in the 
differential system (eqns. BA8B - BA1 II ) are such that cooling 
terms are comparable to terms proportional to oj for funda- 
mental modes (which have the lowest \u\ for a given £), de- 
creasing their relative contribution for higher overtones as w 
increases in magnitude. When combined ir/2 out-of-phase, 
the pressure perturbations from sloshing modes result in two 
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FIG. 3. — Real part of the perturbed azimuthal velocity from linear stability, showing the three-dimensional structure of spiral modes for different combinations 



of £ and m. The top, middle, and bottom rows show Yj +i'Fj , 7 2 + iY 



and y 2 2 + iF 2 " 2 , 



all at t = 0, respectively. From left to right, columns show different 



elevations above the equatorial plane. The fractional stellar size is r, /r s Q = 0.5 for £ = 1 and r, / r s o = 0.6 for £ = 2, so that only the fundamental mode is unstable. 
The equatorial plane of the 1 mode (panel e) has only radial flow. 



almost half-circular regions of opposite polarity oriented to- 
wards the center of the prograde and retrograde fluid streams. 
As the amplitude grows, a protuberance and depression arise 
at opposite sides of the shock, following the pressure pertur- 
bation (©. The overall flow structure rotates with a period 
equal to the oscillation frequency of the sloshing modes. 

From the previous discussion, it can be seen that spiral 
modes need not be limited to the x-y plane, a phase shift of 
7r/2, or to have equal amplitudes. Panel (a) of Figure|2]shows 
the effects of decreasing the phase difference from ?t/2 to 7r/8 
relative to the standard case shown in panel (1) of Figure \T\ 
The spiral pattern is still apparent, though not static anymore. 
Indeed, the animated version of Figure [2] shows that it is in 
fact episodic, arising whenever the two sloshing modes inter- 
act constructively. Panel (b) of Figure [2] shows what happens 
when the amplitudes are not the same: again, the spiral pattern 
is still distinguishable, although more irregular. We take these 
results as an indication that a more general class of spiral- 
like modes can be obtained by combining sloshing modes in 
different ways. This means that a larger region of parame- 
ter space (number of modes, relative phase, relative ampli- 



tude) of initial perturbations can result in spiral-like behavior 
of the shock compared with the limiting case of ±7r/2 phase 
or equal amplitudes. In $4] we show that in fact these modes 
survive in the nonlinear phase, leading to angular momentum 
redistribution below the shock (albeit smaller in magnitude 
than a phase difference of ir/2). 

Adding a third sloshing mode in the z-direction (Y®) results 
in a 45 degree tilt of the "spiral plane" relative to the equa- 
tor. The azimuthal and polar angles determining the normal 
to this plane are set by the phase and amplitude of this mode 
relative to their transverse counterparts, respectively. One can 
thus immediately see that the torque imparted to the star by 
nonlinear spiral modes is in a direction set by the way the in- 
stability is excited. This conclu sion agrees with the results of 
iBlondin & Mezzacappal (120071) . 

In the parameter regime relevant to core-collapse super- 
novae, quadrupolar modes are also significant. As with the 
dipolar case, each of these corresponds to an elementary 
sloshing SASI mode, which can be combined and excited in 
phase to generate a sloshing mode along an arbitrary axis. 
Spiral modes are again generated by exciting at least two of 
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FIG. 4. — Azimuthal velocity in a snapshot of model R5_L1 1P2 at time t = 33fffo. Units and color table are the same as in FigurefJ] and panels from left to right 
correspond to the same altitudes above the equator. The ratio of stellar to shock radius is r*/r s o = 0.5, for which only the fundamental £ = 1 mode is unstable. 
Despite the relatively coarse grid resolution, the agreement with the linear eigenmodes is excellent. The green contour shows the surface with internal pressure 
equal to O.lpiVjL. An animated version of this figure is available in the online version of the article. 

these basic modes out-of -phase. Figure[3]compares the pertur- 
bation to the azimuthal velocity for dipolar and quadrupolar 
spiral modes at different altitudes from the equatorial plane. 
The \m\ = £ modes preserve their spiral structure from pole 
to pole, while the \m\ < £ case has a purely radial flow at 
the equator, and dipolar spiral flows above and below. These 
dipolar spirals are 180 degrees out of phase with each other. 

Increasing the angular degree of the mode results in more 
complicated combinations, although the modes Y^ e always 

have the form sin 6>{cos(^6>),sin(^f?)}. This guarantees the 
existence of two spiral modes with 2£ arms of alternating po- 
larity for all £ > 1. 

4. RESULTS FROM TIME-DEPENDENT SIMULATIONS 

4.1. Linear Growth and Saturation 

In order to trigger an individual spiral mode in our sim- 
ulations, we drop overdense shells with angular dependence 
determined by real spherical harmonics (eqns. l4l- lfTT1l ). We 
locate these shells in the upstream flow in such a way that 
when advected, they arrive at the shock with a time delay that 
corresponds to a relative phase $. In other words, the radial 
positions r\ and ri of a two-shell combination satisfy 




n dr_ 



(12) 



More than two modes are included straightforwardly, as only 
relative phases are relevant. This type of perturbation results 
in the excitation of an isolated spiral mode when the back- 
ground accretion flow does not rotate, allowing comparison 
with eigenfunctions from linear stability analysis, and quan- 
tification of its contribution to the angular momentum redis- 
tribution as it grows exponentially in amplitude. 

Figure 2] shows the azimuthal velocity field at the same al- 
titudes above the equator as in Figure[3]for model R5_L1 1P2, 
for which only the fundamental mode is unstable. Despite 
the coarseness of the grid, the flow structure in the upper row 
of Figure |3] is clearly reproducecfl The online version of the 
article contains an animated version of this figure. 

To test the effects of unstable overtones, we have also 
evolved a configuration with r*/r s o = 0.2, which is closer 
in size to stalled supernova shocks (models R2_Lllf and 
R2_Lllh). In this case 1 =1 has several unstable over - 
tones (see, e.g., Figure 13 of lFernandez & Thompsonl2 009b). 

3 Despite the fact that the amplitudes of the sloshing modes are slightly 
different. 



FIG. 5. — Evolution of the azimuthal velocity in the equatorial plane at 
/ ~ 30fff for models R2_Lllf (left) and R2_Lllh (right), which are closer 
in size to stalled supernova shocks, but have several unstable i = 1 over- 
tones. When the fundamental is excited (left), an irregular spiral is obtained, 
whereas an excitation of the first overtone yields an oblique sloshing mode at 
the fundamental frequency (right). The color scale is the same as in Figuref4] 

which are all excited whenever an I = 1 perturbation is im- 
posed. Figure [5] shows the resulting azimuthal velocity aris- 
ing from perturbations of the type in equation dT2b . employing 
phase differences corresponding to the fundamental and first 
overtones of the Y^ + iY[ l mode. In the first case, shown in the 
left panel, a spiral-like pattern is formed, although the pres- 
ence of unstable overtones makes it irregular. In the second, 
on the right, an initial first overtone spiral becomes overtaken 
by a fundamental sloshing mode in an oblique direction, be- 
cause the latter mode has a larger growth rate. 

As a mode grows in amplitude, the pressure perturbation 
shown in Figure [TJc develops into a protuberance that rotates 
with the oscillation period of the instability. Eventually a 
shock forms at the interface where the two streams move 
towards each other, res ulting in a triple point at the shock 
(Blon din & Shawl 120071) . Figure [6] shows the morphology 
of the shock right before this triple point forms, for models 
R5_L1 1P2, R6_L21P2, and R6_L22P2. 

The initial phase of exponential growth of isolated modes is 
followed by nonlinear mode coupling and saturation. Figures 
[7^ and[7p show the spherical harmonic coefficients 



a, m (t)= / <mY?R s (0,<l>,t), 



(13) 



where R s (8,4>,t) is an isobaric surface tracing the shock, as a 
function of time for models R5_LllP2c and R5_L11P4. In 
both cases, the primary sloshing mode combination remains 
isolated in the linear phase, along with a non-oscillatory in- 
crease in the £ = 2, m = mode due to the increasing oblate- 
ness of the shock, likel y caused by the redistribution of an- 
gular momentum ( 34.21 ). Other modes are excited only after 
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FIG. 6. — Isobaric surfaces (p = 0.1pii'j f0 ), showing the morphology of the shock in the transition to the nonlinear regime for models R5_L1 1P2 at t = 46.5fffo 
(left), R6_L21P2 at ( = 57%o (middle), and R6_L22P2 at t = 52.5fffo (right). In all cases, the ratio r„/r s o is chosen so that only the fundamental mode is unstable, 
for eac h L The meridional line at y = is an artifact of the plotting tool, while the irregularity at the north pole is due to the numerical instability described in 
section 1231 In all cases, the shock deformations grow in amplitude and rotate counterclockwise. An animated version of this Figure is available in the online 
version of the article. 



the primary mode reaches nonlinear amplitudes. Interestingly, 
the first modes to couple are Y^ 2 , which achieve amplitudes 
about 50% lower than the primary mode. The axisymmetric 
£ = 1 mode together with 1 = 2, m = {0,±1} modes couple 
much later, and do not reach as large an amplitude as m = 2. 
We do not show 1 = 3 results to keep Figure [7] legible, but a 
similar trend is observed in that F 3 ±3 couple first and reach the 
largest amplitudes of alH = 3 modes. This sequence of mod e 
coupling differs from that found by Iwaka mi et alj (|2008a), 
who witnessed the 1 = 2, m = 0, ±1 modes to couple first and 
reach the largest amplitude behind I = 1 . However, their non- 
axisymmetric perturbation excites the real spherical harmon- 
ics and Y{ in phase, hence it is not surprising that the mode 
coupling sequence differs. 

To track how the relative phase in the primary spiral mode 
evolves going into the nonlinear phase, we compute a normal- 
ized cross-correlation between the spherical harmonic coeffi- 
cients of the corresponding sloshing modes: 

CC[ai,i,ci_i](A*) = 

/ dt [fli,i(f)-flu] [ai-i(f + AQ-at-i] 

(Jdt J I J dt [ai_i(0-ai,-i] J 

where a^ m is the time average of a^ m and the integrals are 
performed over the same time interval [f m i n ,f m ax]- Figure [8] 
shows the results of applying equation ( TBI to the linear 
and nonlinear phase of modes R5_LllP2c, R5_L11P4, and 
R5_L11P4, which differ only in the initial relative phase of 
the perturbation. In the linear regime, the cross-correlation 
has an initial peak at a time delay Af; = <£>/u; oje . As this 
amounts to shifting the second sloshing mode back into the 
first, equation ( TBI yields values near unitjQ. Subsequent 
peaks arise at At„ = (<fr + mr)/u) osc , (n = 1,2, ...), with an am- 
plitude (— l)"exp(— nTOcigrow/kW); ^ e -' eacn P er i°d down 
in magnitude by the inverse exponential growth rate, with al- 
ternating sign. When equation (fl4l i is applied to the fully sat- 
urated nonlinear phase, the first peak shifts to a larger value of 

4 The fact that it is not exactly unity means that both sloshing modes were 
not excited with the same amplitude, an effect we believe is due to the partial 
dissipation of the second overdense shell when advected by the supersonic 
flow before crossing the shock. 




FIG. 7. — Real spherical harmonic coefficients of the shock (traced by an 
isobaric surface) as a function of time, for models R5_LllP2c (top) and 
R5_L11P4 (bottom). The linear phase lasts approximately until the I = 1 
displacement is ~ 10% of the postshock cavity, or t ~ 30%o- Full saturation 
is reached around t = 60%o ■ Note that the Y^ 1 1 coefficients do not fall into the 
same envelope in the linear phase, indicating that they were not excited with 
the same amplitude. 

Af relative to the linear phase, with no significant change in 
amplitude between subsequent peaks. Furthermore, the three 
nonlinear cross-correlations are nearly identical, with a period 
~ 12?ffo and a time delay ~ 2.5%o, corresponding to a phase 
shift ~ 2tt/5. This suggests that the development of spiral 
modes in the fully saturated stage is independent of the ini- 
tial conditions, for any non-zero initial relative phase between 
sloshing modes. As mentioned in §|3] a larger class of ini- 
tial perturbations can therefore result in spiral modes, not just 
those tuned to a relative phase of tt/2. 
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FIG. 8. — Normalized cross-correlation (eq. 1141 ) applied to models 
R5_LllP2c (top), R5_L11P4 (middle), and R5_L11P8 (bottom). Shown 
is the linear and weakly nonlinear phase (black), covering the time inter- 
val t e [24,50]f ff0 , and the fully saturated phase (red), over t e [100,200]« ff0 
(see FigureQ. The ratio of consecutive peaks in the linear phase is approx- 
imately exp(7Rjg row /aj osc ), or the fraction of an e-folding corresponding to 
one-half of an oscillation period. The initial condition is reflected in the 
cross-correlation taken over the linear phase. However, the time-delay be- 
tween sloshing modes in the nonlinear phase is almost independent of initial 
conditions. 

We have also performed a simulation in which, instead of 
using overdense shells, we impose random cell-to-cell pres- 
sure perturbations below the shock (R5_RAN). For this, we 
employ the setup with r*/r s o = 0.5, for which only the funda- 
mental I = 1 mode is unstable (and the most unstable for all I). 
Figure [9^ shows the corresponding spherical harmonic coef- 
ficients as a function of time. Initially, a sloshing mode along 
the z-axis is triggered due to the form of the grid. However, 
as this mode reaches the nonlinear phase, a sloshing mode 
along the y axis emerges out-of-phase, and then another along 
the x axis follows. Figure |9p shows the normalized cross- 
correlation in the nonlinear phase, for the three different £ = 1 
mode combinations over the time interval t € [200, 300]%), 
showing that all of them are out-of phase. Note that the phase 
delay between Fj and Y[ l is nearly the same as that in Fig- 
ure [8] (modulo half a period), even though is still growing 
in amplitude. T he re sulting angular momentum redistribution 
is discussed in 34.21 

To close this subsection, we address the effects of dimen- 
sionality and resolution on the saturation amplitude. Table [2] 
shows the rms fluctuation of the ai,i spherical harmonic co- 
efficient around its mean for various models. There does not 
appear to be a systematic difference between the saturation 
amplitude in 2D and 3D at this resolution and with this nu- 
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FIG. 9. — Top: spherical harmonic coefficients for model R5_RAN. The 
initial bias towards axisymmetric (Y® and Y®) perturbations is due to the form 
of the grid (more cells closer to the poles). Non-axisymmetric modes do 
indeed grow out of n umerical noise, as expected. Bottom: normalized cross- 
correlation (eq. 1141 ) applied over the time interval t £ [200, 300]%o to the 
I — 1 modes pairwise. All three of them are out-of-phase with each oth er in 
the nonlinear phase, leading to angular momentum redistribution (Figure [T6l . 



TABLE 2 

Saturation Amplitude of y\ mode 



Model Dimensionality Resolution" A«i i b 



r5 LI LR 


2D 


56 x48 


0.505 


r5 LI LRc 


2D 


56 x 48c 


0.424 


r5 LI HR 


2D 


112 x 96 


0.412 


r5 LI HRc 


2D 


112 x96c 


0.454 


R5 Lllx 


3D 


56 x 48 x 96c 


0.473 


R5 L11P2 


3D 


56 x 48 x 96 


0.505 


R5 LllP2c 


3D 


56 x 48 x 96c 


0.499 


R5 L11P4 


3D 


56 x 48 x 96c 


0.504 


R5 L11P8 


3D 


56 x 48 x 96c 


0.481 


R5 Lll HR 


3D 


112x 96x 192c 


0.525 



a The letter c denotes a 5 degree cutout around the polar axis. 
b For models R5_L11P2 and R5_L11_HR the time range em- 
ployed is [100, 147]fflo and [100, 180]f fra , respectively. For all 
others, it is [100,200]%) ■ 



merical method. The changes due to the presence of the cutout 
around the polar axis or a resolution doubling are of the same 
magnitude, ~ 10%. Similarly, whether the sloshing mode is 
isolated or part of a spiral mode seems to make a small differ- 
ence, which is a weak function of the relative ph ase between 
mode s. In the context of parasitic instabilities (Guilet et al. 
2010), this indicates that either (i) the resolution is still too 
low to expose the different growth rates of parasites in 2D and 
3D, or (ii) given the structure of these modes, the difference 
in the growth rates does not translate in sizable differences in 
the saturation amplitudes. 
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4.2. Angular Momentum Redistribution 

Given that the background accretion flow has not net angu- 
lar momentum, any net spin-up of the protoneutron star via 
accreted matter involves the separation of the postshock flow 
into regions with ang ular momentum of opposite sign (e.g., 
iBlondin &" Shaw 2007). In what follows, we explore how the 
linear and nonlinear phases of spiral modes mediate this an- 
gular momentum redistribution, and quantify the magnitude 
of the effect. To ease comparison with other studies, we ex- 
press angular momenta in units of Mr s o 2 , which amounts to 
dividing our resul ts by Ait within the dimensionless unit sys- 
tem described in jj2.ll 

4.2.1. Linear Phase 

Figure [10] shows the z-component of the angular momen- 
tum density integrated over a spherical surface at radius r 



(15) 



h( r ) = r / driprsin6*V0, 



for model R5_L11_HR at various times during the linear 
phase. The upper panel shows times covering a full oscilla- 
tion cycle, and the bottom shows curves rescaled by a fac- 
tor exp [-2w grow (f - fo)] over a longer timescale, with w grow 
the growth rate measured from the spherical harmonic coef- 
ficients. We thus infer that, in the linear phase, (i) a spiral 
mode has a characteristic angular momentum density profile 
which separates the flow into two counter-rotating regions, (ii) 
this profile is non-oscillatory, and (iii) it grows in amplitude 
at nearly twice the growth rate of the mode. The latter is con- 
sistent with the fact that l z (r) is a second-order quantity, as 
the zeroth- and first-order components of equation ( TT3T > van- 
ish. To lowest order in A£/r s o, it is made up of products of 

the form Ap^'Av^ 1 ' and p^Av?, where the superscripts in 
parentheses indicate explicitly the order of the perturbation. 
We can separate out equation $15[ into components 



: /(0)(2) +/ (l)(l) + O(3) 

dflp (0) sin6»Av ( , 2) 



/(0)(2) . 

Z (D(I) = r 3 / dnAp > 



n(!) 



sinflAv^ 



(16) 
(17) 

(18) 



where 0(3) denotes higher order terms. Due to the fact that 
p (0 ' is spherically symmetric, only the I = component of 
sin#Av (2) survives: 



sin6Av ( 2) 


o,o / 







dfiy o °sin0v -O(3). (19) 



The first order components are the real parts of the complex 
perturbations 



sin 0Av« = I £ m [Sv^T^*^ 



(20) 



ov Q£,m L t e 



(21) 



where T'" denotes the usual complex spherical harmonics, 
and the star stands for complex conjugation. Integrating over 




FIG. 10. — Z-component of the ang ular momentum density integrated 
over a spherical surface l z (r) (eq. |15| ) as a function of time for model 
R5_L11_HR. Top: Values at four times covering an entire oi i oscillation 
period, showing that U(r) is non-oscillatory. Bottom: values at a few times 
where the Fj mode is at phase zero, scaled by the square of the mode ampli- 
fication factor relative to tg = 15fffo, exp(-2aj grow — foD, showing that l z (r) 
grows at nearly twice the growth rate of the primary linear spiral mode. 



angles yields 

;(l)(D = 



r 3 ]Tm e llm{ ^< 

Lm 



[lm(Sp Lm )Re(Sva Am ) 



(22) 



-Re(Sp£ !m )Jm(Svn,£,m)] +0(3). 

Figure QT| shows l z (r) at time t = 23.5/fio for model 
R5_L11_HR, where the Tj 1 + iY^ 1 mode is at phase zero (cf. 
Figure |7J. Also shown are 

;(0)(2) 5 with both p (0) and sin0 Ay (2) 

obtained by projection of the corresponding fields fields onto 
Yq, and the £ = 1, m = 1 component of For the latter, 

the real part of the complex amplitude at phase zero is the 
projection onto Y 1 , while the imaginary part is minus the 
component. Most of the angular momentum density comes 
from / <l))(2) , with Z (1)(1) being a ~ 10% correction. Together, 
these terms account for almost all of l,{r), in agreement with 
the fact that modes with £>2,m^Q have a negligible ampli- 
tude in the linear phase. 

For comparison, figures [T2h and b show the Y^ 1 compo- 
nents of the density and azimuthal velocity fields as a func- 
tion of radius at t = 23.5%), for model R5_L11_HR, along 
with the eigenmodes from linear stability analysis normal- 
ized by the value of their respective spherical harmonic co- 
efficients. Reasonable agreement is found between both ap- 
proaches. Thus, most of / <1)(1) can be accounted for by linear 
theory. Figure [T2t shows [sinflvJ^Jco (eq. lO). Note that 
the amplitude is ~ 10 times smaller than the first order £ = 1 
components in Figure [12^. 

The angular momentum redistribution is a smooth process, 
mediated by exponential amplification. Thus, the formation 
of a triple point at the shock is not the main agent behind the 
the spin-up of the inner region, as we have shown that angu- 
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FIG. 1 1. — Angular momentum den sity integrated over a spherical surface 
for model R5_L1 1_HR (red line, eq. 1151 ) as a function of radius, evaluated 
at t = 23.5%o. where the spiral mode is at phase zero. Also shown are the 
dom inant second order components / (0 '< 2 ' and jCwJ (defined in eqns. 1161 - 
fl8l . black and green lines, re spec tively), where the latter is evaluated for the 
t = 1, m = 1 mode only (eq. 1221 ). and their sum (blue dashed lines). The 
vertical dotted line shows the average shock radius. 
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FIG. 12. — Top: I = 1, m = ±1 components of the density at time t = 23.5%o 
for model R5_L11_HR (black and red curves). Also shown are the cor- 
responding values predicted by linear stability theory for a relative phase 
$ = 7r/2, with Fj 1 at zero phase. The amplitude of the eigenfunctions has 
been normalized by the value of the shock spherical harmonic coefficients. 
Middle: same as top, but now for the azimuthal velocity (see Appendix lAlfor 
the definition of Sva). The reasonable agreement with eigenfunctions means 
that /(W) is mostly accounted for by linear theory. Bottom: The projection 
of the azimuthal velocity onto YS, which is second order (eq. 1191 ). Note that 
its amplitude is about 10 times smaller than the first-order I = 1 perturbations 
in panel (b). 
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FIG. 13. — Evolutio n of the three cartesian components of the enclosed an- 
gular momentum (eq. 1231 ) evaluated at r = 0.6r s o, where l z (r) changes sign, 
for £ = 1, m = ±1 spiral modes at different resolution. Shown are L,x, 
and L z for model R5_L1 1_HR (red, black, and blue curves, respectively) and 
R5_LllP2c (green, magenta, and orange curves, respectively). The bulk of 
the angular momentum redistribution occurs during the exponentially grow- 
ing phase, and the saturation value is a weak function of resolution. The 
subsequent evolution shows larger departures. Model R5_L11_HR is only 
evolved until t = 180%o- 

lar momentum redistribution begins as soon as linear modes 
are excited. The triple point may, however, be related to the 
saturation of the primary mode, setting the magnitude of the 
angular momentum redistribution. 

A key element for a predictive theory of the spin-up due 
to exponentially growing spiral modes involves calculation of 
the second order velocity perturbation sinf?v? , since Z (0)(2) 
dominates the spin-up. We surmise that this would r equir e re- 
peating the steps followed for linear perturbations ( lA.ll i. but 
now expanding to second order in A£ /r s o. Combining this re- 
sult wi th a criterion for t he saturation of the SASI (such as that 
from lGuilet et alj|2010l) would allow estimation of the max- 
imum angular momentum redistribution achievable by this 
means. 

4.2.2. Nonlinear Phase and Total Spin-Up 

We focus first on the high resolution I = 1 model, 
R5_L11_HR, and then present results for other parameter 
combinations (Table[T|i. Figure [T3l shows the evolution of the 
three components of the enclosed angular momentum 

h(r) = J h(s)ds (i = x,y,z) (23) 

evaluated at r = 0.6r s o, corresponding to the approximate ra- 
dius at which l z (r) changes sign (see Figure [TOb. Two clear 
phases can be distinguished. First, from t = to 5 Of go, L z 
grows exponentially due to the effects discussed in 94.2.11 
Once other modes start to grow in amplitude, the torque 
becomes oscillatory and saturates together with the primary 
mode. Once full saturation has been reached, the trans- 
verse components L x and L v increase in magnitude and fluc- 
tuate stochastically around zero. The subsequent evolution in 
this fully saturated stage depends somewhat on the resolution 
employed. Model R5_L11_HR seems to remain in a quasi 
steady-state until the end of the integration at t = 180ffto, while 
model R5_L1 lP2c shows large amplitude and long period os- 
cillations in L z , and a secular increase in the magnitude of L x 
and L v at late times. 
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Figure[l4]shows how the angular momentum redistribution 
depends on the polar and azimuthal indices of the primary 
spiral mode, as well as on the relative phase with which it is 
excited. In this case, the ratio r*/r s o is chosen so that only 
the fundamental mode is unstable for either I = 1 or £ = 2, 
hence the primary spiral mode grows isolated in the linear 
phase. The top panel focuses on the low resolution runs that 
excite the £ = 1, m = ±1 modes (R5_L1 lP2c, R5_L1 1P4, and 
R5_L11P8). The maximum spin-up due to the primary spi- 
ral mode (along L z ) is weakly dependent on the initial rela- 
tive phase, and has a magnitude ~ 0.6Mr s o 2 . However, the 
time-delay required to achieve this maximum can differ by a 
factor of two. The transverse components (L x and L y ) do not 
show a significant difference, aside from the secular growth 
in the case $ = 7r/2, which is likely due to resolution effects 
(Figure [TTt. Figure [14b focuses on the £ = 2, m = ±1 modes, 
which still impart a spin-up along the z-axis despite having 
vanishing azimuthal velocity at the equator. Note however 
that the magnitude of the spin-up is much smaller than £ = 1 
spiral modes. The magnitude of the transverse components 
at late times for $ = n/2 are to be taken with caution, again 
due to possible resolution effects. Finally, the £ = 2, m = ±2 
modes show a behavior qualitatively similar to £ = 1, m = ±1, 
although the magnitude of the maximum angular momentum 
redistribution depends strongly on the relative phase. As with 
£ = 2, m = ±1, the maximum spin-up is much smaller in mag- 
nitude than £ = 1 . 

The effects of changing the size of the postshock cavity 
(r*/r s o) on the angular momentum redistribution are shown 
in Figure [15] for model R2_Lllfc. Here, as in Figure^, 
the system has a size more similar to realistic core-collapse 
situations, but is such that many £ = 1 overtones are un- 
stable. Still, a spin-up is observed, which is in fact larger 
than the case r»/r s o = 0.5. The secular growth of the spin- 
up, however, is to be taken with caution due to the low res- 
olution of the model. The angular momentum redistribu- 
tion following the phase of linear growth is still ~ 0.6Mr s o 2 . 
This initial peak, however, comes after a long time-delay of 

130f ff0 ~ 400M[3 /2 (r s0 /150 km) 3 / 2 ms. 

The effect of changing the initial perturbation is shown 
in Figure [16] which shows the angular momentum enclosed 
within r = 0.6r s o for model R5_RAN. The angular momen- 
tum amplification is directly related to the growth of spiral 
modes (compare with Figure[9]l. By the time integration stops, 
Tj 1 is still growing, so L v may increase even more in magni- 
tude. The fact that the spin-up along the z-axis is positive (in 
contrast to models R5_L1 1_HR and R5_L1 lP2c) is a conse- 
quence of the random initial relative phase between the cor- 
responding sloshing modes. The spherical symmetry of the 
problem implies that changing the sign of $ in equation dT2l > 
leads to angular momentum redistribution with the opposite 
sign. 

The order-of-magnitude of the maximum angular momen- 
tum redistribution, Mr s o 2 , can be understood if at saturation, 
the azimuthal velocity i s close to the upstream velocity 
vi . From linear theory (lA.ll i. this follows for shock displace- 
ments of order unity, which is indeed the case for £ = 1 modes. 
Not so straightforward to explain is the size of the fluctua- 
tions in the fully saturated state, which involves interference 
between modes with different amplitudes and phases. The 
basic behavior of the spin-up durin g the exponential growth 
phase agrees with what was seen by Blondin&[Shaw|(2007) 
(their Figure 7), although they did not quantify their results in 
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FIG. 14. — Enclosed angular momentum components at fiducial radial po- 
sitions for models where spiral modes with different spherical harmonic in- 
dices and relative phases are excited. Top: models R5_LllP2c, R5_L11P4, 
and R5_L11P8. The relative phase of the primary spiral mode is related 
to the time required to reach maximum spin-up. Middle: Same for models 
R6_L21P2c, R6_L21P4, and R6_L21P8. This time the angular momentum 
redistribution is more strongly dependent on initial relative phase. The radius 
r = 0.68/" s o is equivalent in depth to that for runs with r*/r s o = 0.5. Bottom: 
Corresponding curves for models R6_L22P2c, R6_L22P4, and R6_L22P8. 
The behavior is somewhat similar to the modes with I = 1, m = ±1, but the 
maximum spin-up does depend on initial phase. Note that the vertical scale 
is not the same. 
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FIG. 15. — Enclosed angular momentum within r = 0.35r s o for model 
R2_Lllfc, for which I = 1 has several unstable overtones (see Figure [5). 
The linear phase is qualitatively similar to that in model R5_L1 lP2c, but the 
maximum angular momentum is larger at longer times. 
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FIG. 16. — Enclosed angular momentum within r = 0.6r s o for model 
R5_RAN (compare with Figure|9j. The appearance of spiral modes induces 
a steady growth in the magnitude of the L x and L z components. Note that the 
model is evolved to later times. 

terms of fundamental physical quanti ties. Our results agree to 
within a factor of a few with those of Blondin & Mezzacappa 
( 120071) . They report typical spin-ups of 2.5 x 10 47 at 250 ms 
(72%) forM= 1.2M and r s0 = 230/1.5 = 153 km, where 1.5 
is our approximate average shock radius in units of r S Q and 
230 km their quoted shock position). At similar times, we find 
0.6Mr s0 2 ~ 10 47 (M/0.36M Q S - 1 )(r s0 /153 km) 2 g cm 2 s" 1 . 
The differences could be caused by several factors: (i) their 
absorbing boundary condition with no cooling, which causes 
the shock to expand with time (e.g.. lBlondin et al1l2003l) . and 
(ii) resolution, which we have already shown to affect the long 
term behavior of the system. 

Assuming that all of the angular momentum enclosed 
within the radius where l z (r) changes sign is accreted onto 
the neutron star, and that the latter has a moment of inertia 
Ins = hs x 10 45 g cm 2 , one can write the minimum period due 
to spiral modes as 



P ~ 8O/45M0 1 , 



150 km 



0.6 

/amp 



T- ms > ( 24 ) 



where / amp is the fraction of Mr s o 2 achieved during the phase 
of exponential growth. If no exponentially growing spiral 
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FIG. 17. — Minimum neutron star spin periods, obtained by assuming 
that all of the angular momentum enclosed within the radius at which l z (r) 
changes sign is accreted onto a star of moment of inertia /ns = 10 45 /45, for 
different models and as a function of time. All components of the angular 
momentum are included. 

mode takes place, one can still expect to achieve / amp < 0.1 
from stochastic fluctuations, in which case the spin period P > 
500 ms. Figure[l7]shows the evolution of the spin period for a 
few representative modes, taking into account the three com- 
ponents of the angular momentum, P = 27t/ns / ^/L 2 + L 2 +Lf, 
with the Lj measured at the radius at which L(r) changes sign. 

This simplistic picture is bound to change when more real- 
istic physics is included. First, because it is not clear yet that 
spiral modes can easily de velop in a convective gain region 
(e.g.. llwakami et al]|2008al) . but also because both the neutri- 
nosphere radius, the mass accretion rate, and the shock radius 
change as a function of time. The onset of explosion further 
complicates the picture, as the mass accretion is cut off and 
the postshock flow expands. Still, anisotropic cold and fast 
downflows piercing the convective region are exp ected to im- 
part s tochastic torques to the forming neutron star (I Tho mpson 
2000), an effect that we cannot capture with our current setup. 

5. SUMMARY AND DISCUSSION 

In this paper we have studied the spiral modes of the 
SASI in the linear and nonlinear regime, when no rotation 
is imposed on the accretion flow, and when neutrino driven 
convection is suppressed. We have combined results from 
linear stability analysis and time-dependent simulations 
to understand the structure and evolution of these non- 
axisymmetric modes. Our main findings follow: 

1. - Spiral modes are most easily understood as two or more 
sloshing modes out of phase. In the linear regime and in the 
absence of rotation, their three-dimensional structure can 
be obtained by linear superposition of known axisymmetric 
eigenfunctions. The parameters that determine a spiral mode 
are the number of sloshing modes involved, their relative 
phases, and amplitudes. Two modes of equal amplitude and 
relative phase equal to ±7r/2 is a limiting case of a more 
general class of non-axisymmetric mode. 

2. - As long as the initial relative phase is not zero, spiral 
modes survive in the nonlinear phase. For the £ = 1, m = ±1 
case, they seem to reach some type of equilibrium relative 
phase once all modes are excited in the fully saturated stage. 
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This equilibrium phase does not seem to depend on the 
initial phase shift. Hence, the range of perturbations needed 
to excite spiral-like behavior leading to protoneutron star 
spin-up is broader than that needed to achieve spiral modes 
with a relative phase of ir/2. The absolute sign of the angular 
momentum component in the inner and outer region depends 
on the sign of the initial relative phase. 

3. - The angular momentum redistribution in the linear 
and weakly nonlinear phase of an isolated spiral mode is 
caused by the spatial dependence of the angular momen- 
tum density, which consists of at least two counterotating 
regions. This division occurs because the eigenmodes have 
at least one radial node in the transverse velocity profile 
(see, e.g., Figure [T2b-c). The angular momentum density is 
non-oscillatory, and increases in magnitude at nearly twice 
the growth rate of the spiral mode. For £ = 1, we have found 
that the dominant term involves a second order perturbation 
to the azimuthal velocity, which couples to the background 
density. We see no apparent relation between the formation 
of the triple point at the shock and the angular momentum 
redistribution, other than a possible role in the saturation of 
the primary mode and hence a cutoff in the growth of the 
spin-up. 

4. - The bulk of the angular momentum redistribution is 
achieved at the end of the phase of exponential growth. 
After 200 dynamical times, our models achieve a maximum 
spin-up of at least 0.6Mr s o 2 along the axis of the primary 
spiral mode, independent of the resolution employed, but 
restricted to £ = I. Modes with I = 2 also lead to spin-up 
along the primary axis, but with a smaller magnitude. 

5. - In the fully saturated stage, where all modes are 
excited due to nonlinear coupling, all components of the 
angular momentum fluctuate with a characteristic magnitude 
< Mr s o 2 / 10. Resolution seems to be relevant for capturing 
the long-term evolution of the spin-up, which can display 
secular increases in magnitude at low resolution. 

Blon din & Mez zacappa ( 2007) found that the development 
of spiral modes at late times is a robust feature of the flow, 
as many different types of perturbation resulted in a simi- 
lar outcome when the accretion flow is non-rotating. Our 
results tend to confi rm this picture . This was not the case 
in the simulations of Iwakami et alj d2009al) . which however 
include neutrino heating and thus convection, fundamentally 
altering the flow dynamics relative to the simpler case with 
no heating. The fact that most progenitor models are likely 
to show some degree of rotation will always make excita- 
tion of a prograde spiral mode more l ikely, as its growth rate 
is larger dYamasaki & Fog lizzo 2008). Indeed, the develop- 
ment of SASI-like spiral modes has been observed in time 
dependent studies of accre tion disks around a Kerr black hole 
( Nag akura & Yamadall2009h 

Are overdense shells with large angular extent realistic? 
Two- and three-dimensional compressible simulations of car- 
bon and oxygen shell burning in a 23M Q star find that density 
fluctuations of order 10% are obtained due to internal waves 
excited in the st ably stratified layers that lie i n between con- 
vective shells dMeakin & Arnettll2006l l2007allbh . These large 
amplitude fluctuations are due to the stro ng stratification, 
and are correlated over large angular scales (Mea krn~& Arnettl 



2007at). The oxygen-carbon interface lies too far out in radius 
to reach the stalled shock during the crucial few hundred mil- 
liseconds after bounce. But if a similar phenomenon were to 
occur above the Si shell, it would have a definite impact on 
the evolution of the stalled shock. Given that typical SASI 
periods are ~ 30 ms, all that is needed is a region of radial 

extent < Vff/cj osc ~ 200M 1 1 ^ 2 r 15 'g 2 km by the time it reaches 
the shock. 

A different question is whether this coherent superposi- 
tion of linear modes can take place in an environment where 
turbulent neutrino-driven convection already operates. The 
convective growth tim e in the gain region is a few ms (e.g., 
iFrver & Youngl l2007) compared to the several tens of ms re- 
quired for the SASI to grow. A critical parameter determining 
the interplay between these two instabilities is the integral of 
the buoyancy fre quency multiplied b y the advection time over 
the gain region (Foglizzo et al. 2006). If this dimensionless 
parameter \ is less than 3, then infinitesimal perturbations do 
not have time to grow before they are advected out of the heat- 
ing region, and a large amplitude perturbation is required to 
trigger convection. Using axisymmetric simulations with ap- 
proximate neutrino transport and a realistic equation of state, 
IScheck et alj ((2008) found that for weak neutrino heating, the 
SASI overturns can actually trigger convection. On the other 
hand, u sing parametric simulations, Fernand ez & Tho mpson 
(2009a) found that when x > 3, convection grows rapidly and 
the vorticity distribution reaches its asymptotic value before 
the SASI achieves significant amplitudes. Two-dimensional 
convection is volume filling, as the vorticity accumulates on 
the largest spatial scales, hence large scale convective modes 
could be exciting dipolar SASI modes. The interplay between 
these two instabilities is not well understood at present, hence 
discussion of spiral modes in this context will have to wait for 
further work. 

A few three-dimensional core-collapse simulations have 
been performed so far, some of these with highly dissipa- 
tive nu merical methods (IFrver & Yo ung 2007; Iwaka mTet alJ 
2008a). These groups find that convection starts on smaller 
angular scales, with convective cells hav ing roughly the size 
of the gain region ( IFrver & Yo ung 2007p, or being mediated 
by hig h-entropy bubbles with a range of sizes ( llwakami et alj 
2008a). In exploding simulations, as the shock expands, con- 
vective cells increase size and a global 1= 1 mode emerges 
(IFrver & Youngll2007t llwakami et aLl l2008a). Persistent spi- 
ral modes do not seem to appear spontaneou sly, they need to 
be explicitly triggered (Iwakami et al. 2008a). 

Very recently, Nordha us et alj (120101) have reported three- 
dimensional, parametric core-collapse simulations using Rie- 
mann solvers and covering the whole sphere. Their results 
are in line with previous studies in that they do not wit- 
ness the development of large scale shock oscillations in the 
stalled ph ase, or spiral modes with no ticeable amplitudes. 
Similarly, Wongwathanarat et al. (2010) have presented full- 
sphere parametric explosions with a Riemann hydrodynamic 
solver, including the contraction of the protoneutron star and 
gr ey neutrino tran s port. Their conclusions are similar to that 
of Nordha us et al.1 (120101) in that no coherent spiral modes are 
observed, with angular momenta saturating at a few times 
10 46 g cm 2 s _I . IFrver & Youngl ((2007) observe specific an- 
gular momenta < 10° cm 2 s -1 imparted stochastically to 
the protoneutron star by anisotropic accretion . These val- 
ues, and those of Wongwatha narat et al] (120101) . are in agree- 
ment with the angular momentum fluctuations we observe in 
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our fully saturated SASI, which on average have a magnitude 
~ 5 x 10 12 (/ amp /0.1)Mo. 3 (r s o/150km) 2 M7.3 cm2 s "'- 
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APPENDIX 

CALCULATION OF LINEAR EIGENMODES 
Background Flow and Perturbation Equations 
The steady accretion flow below the shock is obtained by solving the time-independent Euler equations, using an ideal gas 



equation of state of adiabatic index 7 and the gravity of a point mass M, 

d_ 

dr 

dv r 1 dp 
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The boundary conditions at r = r S Q are given by the Rankine-Hugoniot jump conditions (e.g., lLandau & Lifshitzi ri987). The 
upstream flow is adiabatic and has Mach number 1 at r = r S Q. For a given ratio r* / r s o, the normalization of the cooling function 
is obtained by imposing v r (r*) = 0. 

Foglizzo et al. (2007) introduce Eulerian perturbations of the form shown in equation ([T). To obtain a compact formulation of 
the system, the following perturbation variables are employed 
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corresponding to the perturbed mass flux, energy flux, entropy, and an entropy-vortex combination, respectively. The differential 
system for the radial profile of linear perturb ation s is o btained by p erturbing the time-de pendent fluid equations to first order, and 
rewriting the resulting equations in terms of (IA4b - dA7b , obtaining (Foglizzo et al. 2 007b 
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where /1 2 = 1 -(1 - A4 2 )[£(£+ l)c 2 /r]/u> 2 and A4 is the Mach number. The shock boundary conditions for equations dA8b- (IAl lb 
are obtained from the perturbed shock jump conditions evaluated at the unperturbed shock position r s o (Foglizz o et al.ll2007h . 
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where A£ is the shock displacement, Av = -iwA^ the shock velocity, and the subscripts 1 and 2 refer to quantities above and 
below the shock, respectively. The reflecting boundary condition at r = is 



Sv r ■ 



(l-M 2 ) 



h + SS- 



f 



:0, 



(A16) 



resulting in a complex eigenvalue uj. 



Transverse Velocity Perturbation 
The transverse components of the perturbed Euler equation are 

. . v r dSv r 1 dSp 

-IW6V0+ 7T}T+V r 0W < p + 7T7T=0 
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. v r dSv r 1 
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= 0, 



(A17) 
(A18) 



from which one can readily solve fo r 5vg and Sv^ if the transverse components of the vorticity w are known. Note that, since 
8v r and 5p are proportional to T™ ( j|2.21 i. one has {Sve, 8w$} oc dY™/d9, and {5v^, Swg} oc dY" 1 /(sin9 d(f>). Hence the only 
additional quantity required is the radial profile of the vorticity. 
Taking the curl of the Euler equation yields 



The inertial term can be expanded as 
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where Vq • <5wsi is the angular part of the divergence of the transverse vorticity, <5wq = Swg O+Sw^ <fi. Denoting by B the baroclinic 
term, one obtains 



B r = 

1 

Be = — , 

PQ 



1 dpoddp 1 dpodSp 



B^ = - 



Pa 



rsin9 dr 
1 dpo dSp 
r dr d9 



rsin9 dr 
1 dpo d5p 



r dr d9 



The equation governing the radial profile of the transverse vorticity is then 
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Once the radial profile of <5wq is known, solving for S w r is s traightforward, as no radial d erivati ves ar e invo lved. 

The boundary condition at the shock for equation ( IA24b is obtained from equations ( |A17t and ( IA18K once the transverse 
velocity below the shock is known. Imposing f, • (v2 + <5v2- Vi) = at the shock, where f, are the tangent vectors 
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one obtains (e.g. jFoglizzo et aI1 l2007) 
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where the subscripts 1 and 2 denote values upstream and downstream of the shock, respectively. We then obtain 

5w r . 2 = (A29) 

5w02=-—5v < hT, + —6vr£+——6p2 (A30) 
v r r s0 pr s0 

5w < p2 = —Sv s .2-—Sv r .2 —Sp2, (A31) 

v r r sQ pr s0 

where tildes denote t he rad ia l amp litude (eq. HI). Given that Svg,2 = <^0,2, it follows that Sw^.2 = -Swg.2- With this boundary 
condition, equations ( IA22I ). ( 1A231 I. and ( IA24b imply that Swg(r) = —8w^(f) for all r, and hence 5vg(r) = 5v$(r) as well. For 
convenience, we adopt the notation 

tfvn(r) = 5vg(r) = Sv^(r). (A32) 

Given the radial and angular form of the transverse velocity perturbations, one obtains Sw r (r) = (Vxv) r = everywhere. 

Knowing that t he tr ansverse components of the vorticity have equal and opposite radial amplitudes, one can show from the 
definition of 5K that 



5wg = 



rv r \_l{t+\) p 



(A33) 



that is, th e amp litude of the vorticity can be directly obtained from the basic perturbation variables, without the need for solving 
equation (I A24D . We have nevertheless integrated this equation as a self-consistency check. 

ZEUS-MP IMPLEMENTATION 
Tensor Artificial Viscosity 

Astrophysical hydrodynamic codes that employ finite difference algorithms to solve the Euler equations rely on artificial 
viscosity to broaden shocks over a few cells. This allows differencing of the flow variables without discontinuities, while 
obtaining the correct solutions to the jump conditions away from shocks. The standard practice is to use the prescription of 
IVonNeumann & Richtmverl (Q950), which only takes into account the spatial derivative of the velocity normal to the shock as an 
approximation to the divergence of the velocity field for determining compression. When curvilinear coordinates are used, how- 
ever, this a pproach results in short wavel ength oscillations behind the shock, and spurious heating of homologously contracting 
flows (e.g., Tscharn uter & Winklerll 19791) . To fix this, the artificial viscosity can be formulated as a coo rdinate-invariant tensor, 
which can correctly identify zones where the flow is compressed by becoming active only when V • v < (Tscharnuter & Winkler! 
I1979I) . Additional constraints on the functional form come from requiring that the artificial viscosity doe s not act on homolo- 
gously contracting or shear flows (Tsch arnuter & Winkler 1979; Stone & Norman 1992; Haves et al. 2006). 

In the context of core-c o llapse calculations, a te nsor form of the artificial vi scosity in Zeus-MP has previously been employed 
by llwakami et al.l (l2008al) . Ilwakami et al.l (120 09a). and Iwaka mi et al.1 (l2009bl). They found that unl ess the tensor formulation is 
used, the carbuncle instability (lOuirkll 19941) appears at the poles of the shock (|Iwakami et al.ll2008b | ). 

We have independ ently implemented a tensor artificial viscosity in Zeus-MP along the lines of IStone & Normanl (1 1992b and 
Iwak amTet al.l (l2008al), with some slight mo difications. First, we use the volume difference instead of metric coefficient times 
coordinate difference (IStone & Nor man 1992). This removes singularities at the polar axis whenever the factor sin6d9 appears 
in the denominator. Second, we have made use of the traceless nature of the artificial viscosity tensor to rewrite the ^-component 
of its divergences as 



1 9 < ■ an 233 dsinfl 
rsmOou r sin ft 09 



(V-Q) ( 2) = — -^(^9Q 2 2)--^^r- (Bl) 



(sin 9Q 2 2)+ . Q , (B2) 



r sin 2 6 39 y ' rsintf 36 

where the second expression is that used by IStone & Normanl ( 1 19921) and Ilwakami et al.l d2008al) . This allows to completely 
eliminate the metric coefficients of the denominator, by implementing the volume difference described above. The net effect is 
to reduce the noise in the 9 velocity in the supersonic regions close to the z-axis, without making any significant difference at or 
below the shock. 

When evolving our setup with this artificial viscosity prescription, the flow rem ains spherically symmetric and s mooth for 
an indefinite time as long as there are no perturbations. If, on the other hand, the IVonNeumann & Richtmverl (119501) prescrip- 
tion is used, the carb uncle instability indeed appears out of numerical noise within ~ 10 dynamical times, in agreement with 
Ilwakami et aT](l2008bl) . 

However, when evolving sloshing SASI modes transverse to the z-axis, a runaway sawtooth instability of the <p component of 
the velocity develops around the axis whenever the amplitude (and thus the flow velocity) becomes large. This happens because 
the radial velocity at both sides of the axis is not the same, which in turn is caused by a wake generated by the reflecting axis. The 
reflection of the flow at the axis also compresses it, resulting in enhanced artificial viscosity dissipation at the axis, with irregular 
velocities. A similar phenomenon is observed in modes parallel to the z-axis, although the cause of this is not clear yet. 

The cutout at the polar axis, described in the next subsection, suppresses the growth of this numerical instability to the point 
where meaningful results can be obtained. It does not, however, completely eliminate it. We proceed to the nonlinear phase based 
on the fact that the angular momentum contributed by regions around the shock and close to the polar axis is minor. 
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FIG. 18. — Left: spherical harmonic coefficients for modes Y} and Kr 1 as a function of the half-opening angle 8 mm of the cutout around the polar axis. The 
cases 9 mm = and £> m j n = 5 d eg c orrespond to models R5_L1 1P2 and R5_L1 lP2c, respectively. Right: z-component of the angular momentum density integrated 
over a spherical surface (eq. 1151 ) at t = 24fffo as a function of 6 m m. for the same runs as the left panel. 
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FIG. 19. — Comparison of linear growth rates (left panels) and oscillation frequencies (right panels) from numerical simulations (squares) with values from 
linear stability analysis (dashed lines). The left plot shows results for £ = 1 (r* /r s o = 0.5) and the right i — 2 (r* / r sU = 0.6). Different angular resolutions are color 
coded, with red, green, and blue corresponding to Ng = 48, 96 and 192 cells over the whole range of polar angles, respectively. For each angular resolution, two 
points are shown, one with the baseline radial resolution (such that r s oA9 = Ar at the shock) and another (to the right) with radial resolution doubled. At the 
resolution employed in this study, growth rates are reproduced within 20%, and oscillation frequencies within 10% of the linear stability value. 



Cutout at the Polar Axis 

To investigate the reliability of results obtained with the axis cutout prescription described in § 12.31 we performed two runs 
with half-opening angles 1 and 2.5 degrees, which fill the gap between models R5_L11P2 and R5_LllP2c. Figure [T8k shows 
the spherical harmonic coefficients for modes Yj 1 and Ff 1 as a function of time, for different values of the half-opening angle 
# m i n of the cone that is removed from the grid around the polar axis. Aside from a ~ 10% reduction in the maximum amplitude, 
all curves follow nearly identical trajectories until t ~ 80ffto, where they start to diverge more noticeably. Panel (b) shows 
the angular momentum density integrated over a spherical surface at time t = 24?grj (eq. 1151 . a low -resolution counterpart to 
Figure fTTTi. Again, as the opening angle increases, results decrease by about 10%, consistent with the decrease in amplitude of 
shock oscillations. Other than that, the qualitative behavior is identical. We conclude that our resolution-independent results are 
reliable, with a quantitative uncertainty of at least 10%. 

Linear Growth Rates 

The eigenfrequencies obtained with the method described in Appendix lAl have been verified to high prec ision with time- 
dependent axisymmetric hydrodynamic simulations using the code FLASH2.5 (Fernandez & Thompson 2009b). A basic test for 
the newly implemented setup is thus the reproduction of the linear SASI growth rates from linear stability analysis. This is a 
resolution-dependent requirement, hence convergence tests are needed. Figure[l9]shows growth rates and oscillation frequencies 
as a function of angular and ra dial resolution, for axisymmetric simulations of I = 1 and I = 2 modes. The measurement method 
is the same as that described in Fernand ez & Thompson! (|2009b). With the baseline resolution employed in this study (Ne = 48), 



3D SASI 
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growth rates are captured within ~ 20% and oscillation frequencies within 10%. 
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